<!DOCTYPE html PUBLIC "-//W3C//DTD HTML 4.01 Transitional//EN">
<html><head><meta http-equiv="Content-Type" content="text/html; charset=utf-8"><meta name="generator" content="MATLAB R2016a"><meta http-equiv="X-UA-Compatible" content="IE=edge,IE=9,chrome=1"><title>Immisicble two-phase flow drainage in a subsurface reservoir - example 1 </title><style type="text/css">
* {margin: 0; padding: 0;}
body {text-align: start; line-height: 17.2339992523193px; min-height: 0px; white-space: normal; color: rgb(0, 0, 0); font-family: Consolas, Inconsolata, Menlo, monospace; font-style: normal; font-size: 13.9999990463257px; font-weight: normal; text-decoration: none; white-space: normal; }
h1, h2 {font-weight: normal;}
.content { padding: 30px; }

.S0 { margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S1 { line-height: 26.3999977111816px; min-height: 24px; white-space: pre-wrap; color: rgb(213, 80, 0); font-family: Helvetica, Arial, sans-serif; font-size: 22px; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 15px; margin-right: 10px;  }
.S2 { min-height: 0px; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S3 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }
.S4 { min-height: 0px; color: rgb(0, 95, 206); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S5 { font-family: Helvetica, Arial, sans-serif; margin-left: 0px; margin-top: 10px; margin-bottom: 20px; margin-right: 0px;  }
.S6 { text-align: left; line-height: 20.9999980926514px; white-space: pre-wrap; white-space: pre-wrap; margin-left: 55.9999961853027px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S7 { font-style: italic; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S8 { line-height: 20.576000213623px; min-height: 20px; white-space: pre-wrap; color: rgb(60, 60, 60); font-family: Helvetica, Arial, sans-serif; font-size: 16px; font-weight: bold; white-space: pre-wrap; margin-left: 4px; margin-top: 3px; margin-bottom: 9px; margin-right: 10px;  }
.S9 { margin-left: 3px; margin-top: 10px; margin-bottom: 10px; margin-right: 3px;  }
.S10 { line-height: 15.5926666259766px; min-height: 18px; white-space: nowrap; font-size: 12.6666669845581px; white-space: nowrap; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S11 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(34, 139, 34); white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S12 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S13 { line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 10px; margin-bottom: 9px; margin-right: 10px;  }
.S14 { min-height: 0px; font-family: monospace; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S15 { font-family: monospace; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S16 { min-height: 0px; font-style: italic; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S17 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S18 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(160, 32, 240); white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S19 { color: rgb(64, 64, 64); margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S20 { margin-left: 3px; margin-top: 10px; margin-bottom: 4px; margin-right: 3px;  }
.S21 { text-align: center; line-height: 21px; min-height: 17px; white-space: pre-wrap; font-family: Helvetica, Arial, sans-serif; white-space: pre-wrap; margin-left: 4px; margin-top: 2px; margin-bottom: 9px; margin-right: 10px;  }
.S22 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }
.S23 { line-height: 15.5926675796509px; min-height: 0px; white-space: pre; color: rgb(0, 0, 255); white-space: pre; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 45px;  }
.S24 { min-height: 0px; font-weight: bold; margin-left: 0px; margin-top: 0px; margin-bottom: 0px; margin-right: 0px;  }

.LineNodeBlock {margin: 10px 0 10px 0;}
.LineNodeBlock+.paragraphNode {margin-top: 10px;}
.lineNode {padding-left: 10px; background-color: #F7F7F7; border-left: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.inlineWrapper:first-child .lineNode,.inlineWrapper.outputs+.inlineWrapper .lineNode {padding-top: 5px; border-top: 1px solid #E9E9E9;}
.inlineWrapper:last-child .lineNode,.inlineWrapper.outputs .lineNode {padding-bottom: 5px; border-bottom: 1px solid #E9E9E9;}
.lineNode .textBox {white-space: pre;}
.outputGroup {margin: 2px 0 2px 0; padding: 2px 2px 2px 4px;}
.outputRegion {}
.outputParagraph {color: rgba(64, 64, 64, 1); padding: 10px 0 6px 17px; background: white; overflow-x: hidden;}
.inlineWrapper:last-child .outputParagraph {border-bottom-left-radius: 4px; border-bottom-right-radius: 4px;}
.outputParagraph:empty {margin: 0;}
.inlineElement .symbolicElement {}
.embeddedOutputsSymbolicElement .MathEquation {margin-top: 5px; margin-bottom: 5px;}
.embeddedOutputsSymbolicElement .MathEquation.displaySymbolicElement {margin-left: 15px;}
.embeddedOutputsSymbolicElement .MathEquation.inlineSymbolicElement {margin-left: 5px;}
.embeddedOutputsSymbolicElement {overflow-x: auto; overflow-y: hidden;}
.embeddedOutputsSymbolicElement { overflow: initial !important;}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre; word-wrap: initial; min-height: 18px; max-height: 250px; overflow: auto;}
.textElement {padding-top: 3px;}
.embeddedOutputsTextElement.inlineElement {}
.inlineElement .textElement {}
.embeddedOutputsTextElement,.embeddedOutputsVariableStringElement { max-height: none !important; overflow: initial !important;}
.veSpecifier {}
.veContainer {}
.veSpecifierBox {height: 400px; width: 500px;}
.veSpecifier .veTable {padding-top: 3px; padding-bottom: 4px;}
.veSpecifierBox .veSpecifier .veContainer {position: relative; width: 100%; height: 370px;}
.veSpecifierBox .dijitDialogPaneContent {width: 97% !important; height: 88% !important;}
.veSpecifier .veTable .rowHeadersWrapper {padding-bottom: 0;}
.veSpecifier .veTable .scroller .variableEditorRenderers {padding-right: 3px; -webkit-user-select: none; -moz-user-select: none; -ms-user-select: none;}
.veVariableValueSummary {display: inline-block;}
.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {visibility: hidden; z-index: 0;}
.veMetaSummary {font-style: italic;}
.veSpecifier .veTable .scroller {overflow: hidden;}
.veSpecifier .veTable:hover .scroller {overflow: auto;}
.veSpecifier .veVariableName,.veSpecifier .veDimensionFont {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.veSpecifier .veVariableName {padding-top: 3px;}
.veSpecifier .veDimensionFont {font-style: italic; color: #9A9A9A;}
.veSpecifier .scroller::-webkit-scrollbar-track {background-color: white;}
.veSpecifier .scroller::-webkit-scrollbar-corner {background-color: white;}
.veSpecifier .veTable .topRowHeaderWrapper {border: none; background-color: #F8F9FA;}
.veSpecifier .mw_type_ListBase.showCellBorders,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper,.veSpecifier .veTable .verticalScrollSpacer,.veSpecifier .veTable .horizontalScrollSpacer {border: none;}
.veSpecifier .veTable .dataScrollerNode {border: 1px solid #BFBFBF;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode,.veSpecifier .veTable .dataBody {font-family: Arial; font-size: 13px;}
.veSpecifier .veTable .columnHeaderNode,.veSpecifier .veTable .rowHeaderNode {color: #7F7F7F;}
.veSpecifier .veTable .dataBody {color: #000000;}
.veSpecifier .veTable .columnHeaderNode .cell .drag,.veSpecifier .veTable .columnHeaderNode .cell .header,.veSpecifier .veTable .topHeaderWrapper,.veSpecifier .veTable .bottomRowHeaderWrapper {background-color: #F8F9FA;}
.veSpecifier .veTable .columnHeaderNode .cell .dragBorder {border-right: 1px solid #F8F9FA;}
.veSpecifier .veTable .rowHeaderNode .cell {padding-top: 3px; text-align: center; border-bottom: 1px solid #F8F9FA;}
.veSpecifier .veTable .dataBody .cell .plainText {text-align: right;}
.veSpecifier .veTable .dataBody .row .cell {border-bottom: 1px solid #E9E9E9; border-right: 1px solid #E9E9E9;}
.embeddedOutputsVariableElement {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px; white-space: pre-wrap; word-wrap: break-word; min-height: 18px; max-height: 250px; overflow: auto;}
.variableElement {}
.embeddedOutputsVariableElement.inlineElement {}
.inlineElement .variableElement {}
.variableValue { width: 100% !important; }
.embeddedOutputsMatrixElement {min-height: 18px; box-sizing: border-box; font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.matrixElement .variableValue {white-space: pre; display: inline-block; vertical-align: top; overflow: hidden;}
.embeddedOutputsMatrixElement.inlineElement {}
.embeddedOutputsMatrixElement.inlineElement .topHeaderWrapper {display: none;}
.embeddedOutputsMatrixElement.inlineElement .veTable .body {padding-top: 0 !important; max-height: 100px;}
.embeddedOutputsMatrixElement.inlineElement .veVariableName {padding-top: 0;}
.inlineElement .matrixElement {max-height: 300px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer {white-space: nowrap; padding-top: 10px; margin-bottom: 3px;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis.hide,.embeddedOutputsMatrixElement .matrixElement .verticalEllipsis.hide {display: none;}
.embeddedOutputsMatrixElement .matrixElement .valueContainer .horizontalEllipsis {margin-bottom: -3px;}
.matrixElement { max-height: none !important; }
.dijitTooltipDialog .dijitTooltipContainer .dijitTooltipContents .alertPlugin-alertMessage {min-width: 12px; max-width: 400px; max-height: 300px; overflow: auto;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar {width: 11px; height: 11px;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-track {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-corner {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {border-radius: 7px; background-color: rgba(0, 0, 0, 0.1); border: 2px solid rgba(0, 0, 0, 0); background-clip: padding-box;}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.dijitTooltipDialog .alertPlugin-alertMessage::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb {background-color: rgba(0, 0, 0, 0.1);}
.dijitTooltipDialog .alertPlugin-alertMessage:hover::-webkit-scrollbar-thumb:hover {background-color: rgba(0, 0, 0, 0.2);}
.alertPlugin-alertLine {position: absolute; display: initial; width: 40px; text-align: right; cursor: text;}
.alertPlugin-onTextLine {visibility: hidden;}
.alertPlugin-hasTooltip .alertPlugin-warningImg,.alertPlugin-hasTooltip .alertPlugin-errorImg {cursor: pointer;}
.alertPlugin-isStale {-webkit-filter: opacity(0.4) grayscale(80%); filter: opacity(0.4) grayscale(80%);}
.alertPlugin-alertLine .alertPlugin-errorElement {display: inline-block; margin-right: 4px;}
.alertPlugin-errorImg {position: relative; display: inline-block; background-color: rgb(204, 55, 41); width: 14px; height: 14px; margin-top: 1px; -webkit-border-radius: 7px;}
.alertPlugin-alertLine .alertPlugin-warningElement {display: inline-block; margin-right: 3px;}
.alertPlugin-warningImg {position: relative; display: inline-block; width: 16px; height: 15px; overflow: hidden;}
.diagnosticMessage-wrapper {font-family: Consolas, Inconsolata, Menlo, monospace; font-size: 12px;}
.diagnosticMessage-wrapper.diagnosticMessage-warningType {color: rgb(255,100,0);}
.diagnosticMessage-wrapper.diagnosticMessage-warningType a {color: rgb(255,100,0); text-decoration: underline;}
.diagnosticMessage-wrapper.diagnosticMessage-errorType {color: rgb(230,0,0);}
.diagnosticMessage-wrapper.diagnosticMessage-errorType a {color: rgb(230,0,0); text-decoration: underline;}
.diagnosticMessage-wrapper .diagnosticMessage-messagePart {white-space: pre-wrap;}
.diagnosticMessage-wrapper .diagnosticMessage-stackPart {white-space: pre;}
.embeddedOutputsWarningElement{min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsWarningElement.inlineElement {}
.embeddedOutputsErrorElement {min-height: 18px; max-height: 250px; overflow: auto;}
.embeddedOutputsErrorElement.inlineElement {}
.alertPlugin-warningImg {background-image: url()}
.alertPlugin-errorImg {background-image: url()}
</style></head><body><div class = "content"><div class = 'SectionBlock containment active'><h1 class = "S1"><span class = "S2">Immisicble two-phase flow drainage in a subsurface reservoir - example 1 </span></h1><p class = "S3"><span class = "S2">In this tutorial two-dimensional immiscible two-phase flow drainage of oil by water injection is simulated. Permeability and porosity distributions are extracted from the first layer of the full three-dimensional geological model of </span><a href = 'http://www.spe.org/web/csp/datasets/set02.htm'><span class = "S0">SPE10 comparative solution project</span></a><span class = "S2">. </span></p><p class = "S3"><span class = "S2">This tutorial involves four steps:</span></p><ol class = "S5"><li class = "S6"><span class = "S0">First, permeabiiity and porosity data are loaded from disk files, analysed, and visualised. Next, we verify that these variables are highly correlated as expected from the data.</span></li><li class = "S6"><span class = "S0">Setup auxilliary simulation data, such as wells positions, their flow rates, fluid phases properties, initial pressure and saturation conditions, the time stepping parameters, and the nonlinear saturation solver options.</span></li><li class = "S6"><span class = "S0">Run the time simulation loop involving a sequential implicit solution procedure between pressure and saturation distributions. Next, we analyse and plot the convergence behavior and the performance of the nonlinear saturation solver.</span></li><li class = "S6"><span class = "S0">Post-process of the numerical results by plotting the spatial distributions of simulated pressure and saturation at the simulation end. Different presentation graphics styles such as </span><span class = "S7">mapping</span><span class = "S0"> and </span><span class = "S7">contouring</span><span class = "S0"> modes are explored. </span></li></ol></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S8"><span class = "S2">Load and plot the petrophysical (permeability and porosity) data</span></h2><p class = "S3"><span class = "S2">We considere a simple one layer rectangular domain. In this case we explicitly set the domain lengths along each space dimension as follows</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% Domain size along X, Y &amp; Z directions</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Dx = 365.76; Dy = 670.56; Dz = 0.6096;</span></p></div></div><p class = "S13"><span class = "S2">Next, we set the number of cells along each space direction and the total number of grid cells. Next, we add those variables into the </span><span class = "S14">Grid</span><span class = "S2"> structure</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Nx = 60;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Ny = 220;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Nz = 1; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">N  = Nx*Ny*Nz;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.Nx = Nx; Grid.Ny = Ny; Grid.Nz = Nz; Grid.N = N;</span></p></div></div><p class = "S13"><span class = "S2">To fullfill the need of internal compuattional tasks we add the uniform grid spacings, </span><span class = "S14">hx</span><span class = "S2">, </span><span class = "S14">hy</span><span class = "S2">, and </span><span class = "S14">hz</span><span class = "S2">, along the respective directions to the </span><span class = "S14">Grid</span><span class = "S2"> structure</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.hx = (Dx/Nx); Grid.hy = (Dy/Ny); Grid.hz = (Dz/Nz);</span></p></div></div><p class = "S13"><span class = "S2">Additionally, we need to setup other few variables for the </span><span class = "S14">Grid</span><span class = "S2"> structure such as, </span><span class = "S14">V</span><span class = "S2"> (the uniform cells volume), and </span><span class = "S14">compr</span><span class = "S2"> (the total compressibility). Notice that the compressibility data is a 3D shape array whose size conforms to the grid dimensions. The MATLAB command </span><a href = 'https://fr.mathworks.com/help/matlab/ref/ones.html'><span class = "S15">ones</span></a><span class = "S14">(Nx,Ny,Nz)</span><span class = "S2"> constructs such an array with a unit values for all cells, and thus by multiplying it with a constant values we obtain a uniform distribution of the compressibility data over the computational domain. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.V     = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.compr = 4.4e-10.*ones(Nx,Ny,Nz);</span></p></div></div><p class = "S13"><span class = "S2">Now, we load porosity data from a previously saved mat file '</span><span class = "S16">spe_layer1_phi.mat</span><span class = "S2">' which is included in the data subfolder of the included files. porosity data is stored in </span><span class = "S16">phi_layer</span><span class = "S2"> variable (because we can store more than one array in a mat file).</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">load </span><span class = "S18">'spe_layer1_phi.mat' phi_layer</span><span class = "S12">;</span></p></div></div><p class = "S13"><span class = "S2">Now, we will check if the loaded data is within acceptable bounds. A situation which occur many often in practice. First, look if there exists porosity values below or close to 0+. This can be simply done with the following command:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S12">length(find(phi_layer&lt;eps))</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsVariableElement" style="width: 658px;"><div class="variableElement">ans = 68</div></div></div></div></div><p class = "S13"><span class = "S2">indicating that the porosity is below the smallest floating point precision number </span><a href = 'https://fr.mathworks.com/help/matlab/ref/eps.html'><span class = "S15">eps</span></a><span class = "S2"> (effective zero on the target computer) at 68 cells. To correct for this misleading data we set these values to a prescribed minimum porosity threshold (1e-3 herein). Otherwise, the saturation solver will fail since dividing by zero porosity will result in an overflow.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">phi_layer = max (phi_layer(:), 1e-3);</span></p></div></div><p class = "S13"><span class = "S2">The </span><span class = "S14">Grid</span><span class = "S2"> structure is promoted with a new variable </span><span class = "S14">por</span><span class = "S2"> conforming to the </span><span class = "S14">Grid</span><span class = "S2"> shape. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.por = reshape(phi_layer',Nx,Ny,Nz); </span></p></div></div><p class = "S13"><span class = "S2">Finally, permeability data is loaded from the mat file '</span><span class = "S16">spe_layer1_perm.mat</span><span class = "S2">' and variable </span><span class = "S16">perm_layer</span><span class = "S2">. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">load </span><span class = "S18">'spe_layer1_perm.mat' perm_layer</span><span class = "S12">;</span></p></div></div><p class = "S13"><span class = "S2">Tha data units is converted into SI unit of the permeability (m^2) because it is given in Dary units in the SPE10 dataset. Again, the </span><span class = "S14">Grid</span><span class = "S2"> structure is promoted with a new variable </span><span class = "S14">K</span><span class = "S2"> conforming to the </span><span class = "S14">Grid</span><span class = "S2"> shape but with 3 components since the permeability data assumes an anisotropic diagonal tensor. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.K = perm_layer*9.869232667160130e-13;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">Grid.K = reshape(Grid.K',3,Nx,Ny,Nz); </span></p></div></div><p class = "S13"><span class = "S2">Now, we're ready the plot the porosity and permeability data to see how they look like. First, we will plot the porosity data in linear scale. Next, we plot permeability along X direction in log10 logarithmic scale to depict its extent over several orders of magnitude. </span></p><p class = "S3"><span class = "S2">To do this, a cartesian grid </span><span class = "S14">G</span><span class = "S2"> data structure is setup from the domain extent ans spacings. This will help to do the plotting commands more efficiently.   </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);</span></p></div></div><p class = "S13"><span class = "S2">Start a new MATLAB figure and specify the colormap corrsponding the color palette to use and the number of colors in this palette.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">figure;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">colormap(jet(16));</span></p></div></div><p class = "S13"><span class = "S2">Plot the porosity data into a </span><a href = 'https://fr.mathworks.com/help/matlab/ref/subplot.html'><span class = "S15">subplot</span></a><span class = "S2"> in the same window but at le left hand side. The </span><span class = "S14">PlotCellData</span><span class = "S2"> function plots a map of the variable </span><span class = "S14">Grid.por</span><span class = "S2"> using information in the structure variable </span><span class = "S14">G</span><span class = "S2">. The '</span><span class = "S14">EdgeColor</span><span class = "S2">' option is set to '</span><span class = "S14">none</span><span class = "S2">' to hide the display o the gridlines of the underlying grid. Since the grid is so dense in our modl, letting this option active the resulting plot will be obsured. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,1);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">PlotCellData(G, Grid.por(:), </span><span class = "S18">'EdgeColor'</span><span class = "S17">,</span><span class = "S18">'none'</span><span class = "S12">);</span></p></div></div><p class = "S13"><span class = "S2">next, this plot x and y axis are set as dependent and we add a color bar and a title. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">axis </span><span class = "S18">tight equal</span><span class = "S17">; colorbar, title(</span><span class = "S18">'Porosity map [-]'</span><span class = "S12">);</span></p></div></div><p class = "S13"><span class = "S2">The same sequence of commands are repeated for the right plot to obtain a map of Kx in log10 scale. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,2);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">PlotCellData(G, reshape(log10(Grid.K(1,:,:,:)),N,1), </span><span class = "S18">'EdgeColor'</span><span class = "S17">,</span><span class = "S18">'none'</span><span class = "S12">);</span></p></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S17">axis </span><span class = "S18">tight equal</span><span class = "S17">; colorbar, title(</span><span class = "S18">'log_{10}(K_x) permeability map [m^2]'</span><span class = "S12">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 658px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S13"><span class = "S2">Clear phi_layer variable from memory space since it's not needed further. This is more relevant for large scale applications or on computers with very modest ressources. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">clear </span><span class = "S18">phi_layer</span><span class = "S12">;</span></p></div></div><p class = "S13"><span class = "S2">Now, we check if the log10-permeability and porosity fields ara correlated somehow. First, let's start with the first component of the diagonal permability tensor, Kx:</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S12">corrcoef(reshape(log10(Grid.K(1,:,:,:)),N,1),reshape(Grid.por(:,:,:),N,1))</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsMatrixElement" data-width="658" style="width: 658px;"><div class="matrixElement veSpecifier"><div class="veVariableName"><span>ans = </span><span class="veVariableValueSummary veMetaSummary">2x2 double</span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72.5,&quot;totalColumns&quot;:2,&quot;totalRows&quot;:2,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 147px;">    1.0000    0.8819
    0.8819    1.0000
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div><p class = "S13"><span class = "S2">Off-diagonals of the correlation matrix are high positive values indicating a good positive correlation between the two variables. Let's now do the same thing for the y and z components of the permeability to conclude on their correlation with the porosity. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S12">corrcoef(reshape(log10(Grid.K(2,:,:,:)),N,1),reshape(Grid.por(:,:,:),N,1))</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsMatrixElement" data-width="658" style="width: 658px;"><div class="matrixElement veSpecifier"><div class="veVariableName"><span>ans = </span><span class="veVariableValueSummary veMetaSummary">2x2 double</span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72.5,&quot;totalColumns&quot;:2,&quot;totalRows&quot;:2,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 147px;">    1.0000    0.8819
    0.8819    1.0000
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S12">corrcoef(reshape(log10(Grid.K(3,:,:,:)),N,1),reshape(Grid.por(:,:,:),N,1))</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsMatrixElement" data-width="658" style="width: 658px;"><div class="matrixElement veSpecifier"><div class="veVariableName"><span>ans = </span><span class="veVariableValueSummary veMetaSummary">2x2 double</span></div><div class="valueContainer" data-layout="{&quot;columnWidth&quot;:72.5,&quot;totalColumns&quot;:2,&quot;totalRows&quot;:2,&quot;charsPerColumn&quot;:10}"><div class="variableValue" style="width: 147px;">    1.0000    0.7651
    0.7651    1.0000
</div><div class="horizontalEllipsis hide"></div><div class="verticalEllipsis hide"></div></div></div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S8"><span class = "S2">Setup auxilliary simulation data</span></h2><p class = "S3"><span class = "S2">We further add a five-spot well pattern composed from 1 injector and 4 producers. Production wells are placed at the domain corners while the injection well is located at the center. A quarter of the volumetric injection flow rate, </span><span class = "S14">Qi</span><span class = "S2">, is produced at each of the four production wells. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% define some useful time units conversion constants </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">second = 1; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">minute = 60*second;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">hour   = 60*minute; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">day    = 24*hour; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">week   = 7*day;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">month  = 30*day; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">year   = 365*day;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12"></span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% Injection/production flow rates</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">Qi         = 25/hour;                  </span><span class = "S11">% total flow rate (converted from m^3/h to m^3/s)</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">prod_wells = [1,Nx,N-Nx+1,N];          </span><span class = "S11">% grid indices of production wells</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">Qw         = zeros(N,1);               </span><span class = "S11">% initialize production/injection flow rates array</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">Qw(prod_wells) = -Qi/4*ones(1,4);      </span><span class = "S11">% flow rates at production wells </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">inj_well   = (Nx*(Ny-1))/2;            </span><span class = "S11">% grid index of injection well</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">Qw(inj_well)   = +Qi;                  </span><span class = "S11">% injection flow rate</span></p></div></div><p class = "S13"><span class = "S2">Notice that (by convention), flow rates have to be positive for injection wells and negative for production wells. </span><span class = "S14">Qw</span><span class = "S2"> is a distributed 1xN array which holds distributed flow rates. It has zero values everywhere except at injection/production cells. </span></p><p class = "S3"><span class = "S2">Next, we initialize two fluid objects </span><span class = "S14">oil</span><span class = "S2"> and </span><span class = "S14">brine</span><span class = "S2"> in which fluid properties data corresponding to a more dense and viscous </span><span class = "S16">Oil</span><span class = "S2"> phase and less viscous and dense </span><span class = "S16">Saltwater</span><span class = "S2"> phase, respectively. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% fluid properties units: density [Kg/m^3], viscosity [Pa.s], residual saturation []</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">oil   = Fluid(</span><span class = "S18">'Oil'</span><span class = "S17">,      [1200, 3e-3,  0.3]);   </span><span class = "S11">% resident fluid</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">brine = Fluid(</span><span class = "S18">'Saltwater'</span><span class = "S17">,[1030, 5e-4, 1e-2]);   </span><span class = "S11">% injected fluid</span></p></div></div><p class = "S13"><span class = "S2">Initialize fluid global pressure in the computational domain to a constant value of 100 bars and the initial brine saturation to its residual value (i.e. </span><span class = "S14">brine.Sr</span><span class = "S2">)</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% define some useful pressure units conversion constants</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">Pascal = 1;                                      </span><span class = "S11">% Pascal SI unit </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">bar    = 1e+5*Pascal;                            </span><span class = "S11">% bar</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">psi    = 7.042253521126761e+03*Pascal;           </span><span class = "S11">% UK pressure unit = lbf/in^2</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">atm    = 1.013249965828145e+05*Pascal;           </span><span class = "S11">% atmosphere eqv unit </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12"></span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">P = 100*bar*ones(Nx,Ny,Nz);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">S = brine.Sr.*ones(N,1);</span></p></div></div><p class = "S13"><span class = "S2">Define time stepping parameters including total simulation time, </span><span class = "S14">T</span><span class = "S2">, the number of time steps, </span><span class = "S14">nt</span><span class = "S2">, and the global time step, </span><span class = "S14">dt.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">T   = 50*day;                      </span><span class = "S11">% 50 days</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">nt  = 25;                          </span><span class = "S11">% number of time steps</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">dt  = T/nt;                        </span><span class = "S11">% time step</span></p></div></div><p class = "S13"><span class = "S2">For the nonlinear solvers we need to pass an optional structure with options holding the absolute tolerance, tol, and the maximum number of allowed Newton-Raphson iterations, maxiter. The latter corresponds to the inner iteration count since the most outer loop corresponds to time stepping. The last option, min_dt, is the minimum allowed time step size for the saturation solver. This is needed because, as explained later in more detail, we use an adaptive time stepping algorithm which cancels out if the time step size becomes smaller than this prescribed threshold. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">opt.tol     = 1e-5;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">opt.maxiter = 100;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">opt.min_dt  = 30*second;</span></p></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S8"><span class = "S2">Run the time simulation loop - Sequential implicit solution</span></h2><p class = "S3"><span class = "S2">Now, we will build the time loop involving a sequential solution of the governing global pressure equation and the nonlinear non-wetting phase saturation equation formulated in fractional form. The pressure and saturation solvers are both implicit implying an unconditional stability regardlless of the magnitude of the selected time step. </span></p><p class = "S3"><span class = "S2">The pressure equation, neglecting the gravitational and capillary pressure terms is given as:</span></p><p class = "S21"><span style="vertical-align:-15"><img src="" width="182.5" height="34.5" /></span><span class = "S2"> </span></p><p class = "S3"><span class = "S2">where </span><span style="vertical-align:-7"><img src="" width="11.5" height="22" /></span><span class = "S2"> is the total compressibility, </span><span style="vertical-align:-3"><img src="" width="12.5" height="18" /></span><span class = "S2"> is the rock porosity, </span><span style="vertical-align:-3"><img src="" width="11" height="18" /></span><span class = "S2"> is the fluid pressure, </span><span style="vertical-align:-7"><img src="" width="75.5" height="22" /></span><span class = "S2"> is the total mobility given as a sum of oil and water mobilities, respectively. A fluid phase mobility </span><span style="vertical-align:-7"><img src="" width="12.5" height="22" /></span><span class = "S2"> is defined as the ratio of its relative permeability </span><span style="vertical-align:-7"><img src="" width="18" height="22" /></span><span class = "S2"> to its viscosity </span><span style="vertical-align:-7"><img src="" width="13.5" height="22" /></span><span class = "S2">. </span><span style="vertical-align:-7"><img src="" width="16" height="22" /></span><span class = "S2"> is the water injection/withdawal mass flow rate [M/T] and </span><span style="vertical-align:-7"><img src="" width="18" height="22" /></span><span class = "S2"> is the density of the water phase [M/L^3]. Threfore </span><span style="vertical-align:-15"><img src="" width="15.5" height="34.5" /></span><span class = "S2"> is the volumetric flow rate.</span></p><p class = "S3"><span class = "S2">The saturation equation is a nonlinear diffusion-convection like equation for the water phase given as:</span></p><p class = "S21"><span style="vertical-align:-15"><img src="" width="134.5" height="34.5" /></span></p><p class = "S3"><span class = "S2">where </span><span style="vertical-align:-7"><img src="" width="15.5" height="22" /></span><span class = "S2"> is the water phase saturation, </span><span style="vertical-align:-7"><img src="" width="36.5" height="22" /></span><span class = "S2"> is the saturation-dependent fractional flow function for the water phase, </span><span style="vertical-align:-7"><img src="" width="76.5" height="22" /></span><span class = "S2"> is the total velocity (i.e. the sum of fluid velocities of each phase). </span></p><p class = "S3"><span class = "S2">Herein, the pressure solver is linear but requires successive updating of the total mobility of both fluid phases which is the task of </span><span class = "S14">TwoPhasePressure </span><span class = "S2">function. The latter wraps the </span><span class = "S16">Two-Point Flux Approximation solver</span><span class = "S2"> (</span><span class = "S14">TPFA</span><span class = "S2"> function) such that </span><span style="vertical-align:-7"><img src="" width="38.5" height="22" /></span><span class = "S2"> is iteratively updated from spatial distribution of water saturation calculated in the previous time step. </span></p><p class = "S3"><span class = "S2">The saturation solver uses an heuristic adaptive time steeping procedure combined with an inner Newton-Raphson (NR) iteration loop. During the inner NR iterations, the sparse system Jacobian matrix </span><span style="vertical-align:-15"><img src="" width="67" height="34.5" /></span><span class = "S2"> is evaluated analytically to accelerate the algorithm convergence. If the NR iteration fails the current time step is subdivided (i.e. the number of local time steps is multiplied by 2) and the whole time step is repeated again from the beginning. The latter mechanism is called back-stepping. The local NR failure is due to either (i) an insuffisant or small number of allowed maximal iterations to achieve the required accuracy. Hence, when setting the NR options the user should keep a good balance between the </span><span class = "S14">maxiter</span><span class = "S2"> and </span><span class = "S14">tol</span><span class = "S2"> options; (ii) second, to a high or very difficult level of prescribed accuracy which cannot be achived for badly preconditioned problems. So, the absolute tolerance should be kept to reasonable values (in practice values between 1e-3 and 1e-4 are suggested); (iii) and finally, when the selected global time step results in large local pressure jumps, the same time step is generally inappropriate for the saturation solver implying a smaller time step advancement of the saturation solver to absorp the sudden numerical shock caused by the large pressure jump. </span></p><p class = "S3"><span class = "S2">In general, this adaptive time stepping strategy is effective in practice. Very often the algorithm starts with smaller saturation time steps at the begining of the simulation when initial conditions are far from the expected solution. This sub-time step will gradually increase in the next time levels if there is no abrupt changes in the boundary conditions (i.e. wells pumping or production schedules) or materials and fluids properties. </span></p><p class = "S3"><span class = "S2">Let's run this two-phase flow oil drainage simulation</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">tic;              </span><span class = "S11">% start timing the time marching loop</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S22">for </span><span class = "S12">t=1:nt</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   fprintf(</span><span class = "S18">'\n'</span><span class = "S12">); </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   fprintf(</span><span class = "S18">'Solving two-phase flow problem. Time = %f days\n'</span><span class = "S12">, t*dt/day);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">   </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   </span><span class = "S11">% call TPFA flow solver adapted to two-phase flow </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">   [P,V] = TwoPhasePressure(Grid,S,brine,oil,Qw,P,dt); </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">   </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   </span><span class = "S11">% solve for brine saturation </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   [S,</span><span class = "S17">out</span><span class = "S12">(t)] = ImplicitSaturation(Grid,S,brine,oil,V,Qw,dt,opt);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">   </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   </span><span class = "S11">% quit any time step if the saturation solver fails </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">   </span><span class = "S22">if </span><span class = "S17">out(t).convergence==false, </span><span class = "S22">break</span><span class = "S17">; </span><span class = "S23">end</span></p></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S23">end</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 2.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">...........<br>Converged in 11 time backsteps and 4285 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 4.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">....<br>Converged in 4 time backsteps and 82 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement">Solving two-phase flow problem. Time = 6.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">....<br>Converged in 4 time backsteps and 94 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 8.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">....<br>Converged in 4 time backsteps and 97 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 10.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">..<br>Converged in 2 time backsteps and 35 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 12.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">..<br>Converged in 2 time backsteps and 29 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 14.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">..<br>Converged in 2 time backsteps and 29 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 16.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">.<br>Converged in 1 time backsteps and 14 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 18.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">.<br>Converged in 1 time backsteps and 14 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 20.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">.<br>Converged in 1 time backsteps and 15 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement">Solving two-phase flow problem. Time = 22.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 9 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 24.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 7 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 26.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 5 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 28.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 7 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 30.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 5 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 32.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 6 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 34.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 10 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement">Solving two-phase flow problem. Time = 36.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 7 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 38.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 6 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 40.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 6 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 42.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 6 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 44.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 7 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 46.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 6 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 48.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Converged in 0 time backsteps and 6 Newton-Raphson iterations</div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement"></div></div><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Solving two-phase flow problem. Time = 50.000000 days</div></div><div class="inlineElement embeddedOutputsTextElement scrollableOutput" data-width="658" style="width: 658px; max-height: 29px;"><div class="textElement">Converged in 0 time backsteps and 5 Newton-Raphson iterations</div></div></div></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S12">toc;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsTextElement" data-width="658" style="width: 658px;"><div class="textElement">Elapsed time is 23.195861 seconds.</div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S2">From the convergence output we notice that this sequential solution algorithm went very smoothly. There was a difficulty to reach convergence in the first time step because it seems that 2 days time step was very large for the saturation solver. Hence, the initial time step was automatically subdivided 11 times resulting to an effective saturation time step close to 84 seconds as shown below from the fields of the returned output structure </span><span class = "S14">out</span><span class = "S2"> for the first time period. The </span><span class = "S14">dt</span><span class = "S2"> field returns the selected time step by the saturation solver during this period. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S12">out(1)</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsVariableStringElement" data-width="658" style="width: 658px;"><div class="textElement"><div>ans = </div><div>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;convergence: 1<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;backsteps: 11<br>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp; dt: 84.375000000000000<br>&nbsp;&nbsp;&nbsp;&nbsp;nr_iterations: 4285<br><br></div></div></div></div></div></div><p class = "S13"><span class = "S2">Now, we plot a summary of the saturation solver performance during the simulation. In the left side of the figure we plot the evolution of the internally selected time step versus time, while in the right figure the evolution of the total number of inner Newton-Raphson iterations versus time is plotted. Both figure are best depicted using a </span><a href = 'https://fr.mathworks.com/help/matlab/ref/semilogy.html'><span class = "S15">semilogy</span></a><span class = "S2"> axis. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12"></span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,1); </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S22">for </span><span class = "S17">t=1:nt, </span><span class = "S17">dt_history</span><span class = "S17">(t) = out(t).dt; </span><span class = "S23">end </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">semilogy((1:nt)*dt/day, dt_history/day, </span><span class = "S18">'-ro'</span><span class = "S12">); </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">title(</span><span class = "S18">'Evolution of the saturation time-step'</span><span class = "S17">); xlabel(</span><span class = "S18">'Time [days]'</span><span class = "S17">); ylabel(</span><span class = "S18">'Time step [days]'</span><span class = "S17">); grid </span><span class = "S18">on</span><span class = "S12">;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12"></span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,2); </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S22">for </span><span class = "S17">t=1:nt, </span><span class = "S17">nr_history</span><span class = "S17">(t) = out(t).nr_iterations; </span><span class = "S23">end </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">semilogy((1:nt)*dt/day, nr_history, </span><span class = "S18">'-b*'</span><span class = "S12">); </span></p></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S17">title(</span><span class = "S18">'Newton-Raphson iterations'</span><span class = "S17">); xlabel(</span><span class = "S18">'Time [days]'</span><span class = "S17">); ylabel(</span><span class = "S18">'Total NR iterations'</span><span class = "S17">); grid </span><span class = "S18">on</span><span class = "S12">;</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 658px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S13"><span class = "S2">We observe thet the internal saturation time step gradually increases to 1.5 hour, 6 hours, 12 hours and finally to 2 days starting from 22 days after the simulation start. So, afterwards the pressure and saturation time steps becomes equal meaning that the saturation extent of the plume becomes significant enough to amortize the pressure jumps immediately. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S24">Post-processing of the numerical solution - mapping mode</span></p><p class = "S3"><span class = "S2">Since the main simulated variables are pressure </span><span style="vertical-align:-3"><img src="" width="11" height="18" /></span><span class = "S2"> and water saturation </span><span style="vertical-align:-7"><img src="" width="15.5" height="22" /></span><span class = "S2"> we will plot them directly using the supplied functions which make use of MATLAB's excellent built-in graphics capabilities without writing the results into intermediate text or binary files in a required format and loading them back from an external third-party post-processing software. So, with very few commands we get directly the results ready for quality publication freeing the modeller from all intermediate tasks and hence speeding the data-processing-simulation-post-processing workflow significantly. </span></p><p class = "S3"><span class = "S2">We will use here the routine </span><span class = "S14">PlotCellData</span><span class = "S2"> to produce the figures in </span><span class = "S16">mapping</span><span class = "S2"> mode where each cell in the grid is plotted with a constant color corresponding to a linearly interpolated value between the upper and lower colors in the used color palette. This palette fits the range between the maximal and minimal values detected in the data array.  </span></p><p class = "S3"><span class = "S2">The </span><a href = 'https://fr.mathworks.com/help/matlab/ref/figure.html'><span class = "S15">figure</span></a><span class = "S2"> command without any argument opens a new MATLAB figure (or axis in MATLAB's terminology) so that the previously opened figures are preserved and all following drawing commands will be placed in this new figure. The </span><a href = 'https://fr.mathworks.com/help/matlab/ref/colormap.html'><span class = "S15">colormap</span></a><span class = "S2"> command is optionally called before the drawing commands to choose among the available color palettes, here we choose the </span><a href = 'https://fr.mathworks.com/help/matlab/ref/jet.html'><span class = "S15">jet</span></a><span class = "S2"> palette with 256 colors distributed between blue and red. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">figure; </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">colormap(jet(256));</span></p></div></div><p class = "S13"><span class = "S2">Next, we plot the pressure after conversion from </span><span class = "S14">Pascal</span><span class = "S2"> to </span><span class = "S14">bar</span><span class = "S2"> units. We call the </span><a href = 'https://fr.mathworks.com/help/matlab/ref/subplot.html'><span class = "S15">subplot</span></a><span class = "S2"> command to specify that we want to make a plot matrix with 1 row and 2 columns in the target figure. The last argument in the subplot command is the active subplot to which the next plotting commands are directed. The </span><span class = "S14">PlotCellData</span><span class = "S2"> function takes additional property-value pairs (à la MATLAB). These options are exactly those supported by the low level </span><a href = 'https://fr.mathworks.com/help/matlab/ref/patch.html'><span class = "S15">patch</span></a><span class = "S2"> command to control the internal details of the plot. To study all these options type </span><span class = "S14">help patch</span><span class = "S2"> command or go to </span><a href = 'https://fr.mathworks.com/help/matlab/ref/patch.html'><span class = "S15">patch</span></a><span class = "S2"> window in the online help system.  Here the '</span><span class = "S14">EdgeColor</span><span class = "S2">' property specifies the color of the underlying grid edges, by setting this property to '</span><span class = "S14">none</span><span class = "S2">' they will be hidden in the current plot. In this example the grid is very dense. So, drawing all grid edges will obscure the reulting plots but in other applications where the grid is coarse this might not be necessary.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% plot fluid pressure on left side of current figure</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,1);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">PlotCellData(G, P(:)/bar, </span><span class = "S18">'EdgeColor'</span><span class = "S17">,</span><span class = "S18">'none'</span><span class = "S12">);</span></p></div></div><p class = "S13"><span class = "S2">By default MATLAB plots fit to the rectangular window inside which they're displayed. A good idea here is to call the </span><a href = 'https://fr.mathworks.com/help/matlab/ref/axis.html'><span class = "S15">axis</span></a><span class = "S2"> command to make x and y axis scales dependent using </span><span class = "S14">equal </span><span class = "S2">option. Likewise, </span><span class = "S14">tight</span><span class = "S2"> wil readjust the plot to fit at best to actual window after axis scaling. Finally, we decorate the plot with a </span><a href = 'https://fr.mathworks.com/help/matlab/ref/colorbar.html'><span class = "S15">colorbar</span></a><span class = "S2"> and a custom </span><a href = 'https://fr.mathworks.com/help/matlab/ref/title.html'><span class = "S15">title</span></a><span class = "S2">. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">axis </span><span class = "S18">tight equal</span><span class = "S17">; colorbar; title(</span><span class = "S18">'Pressure [bar]'</span><span class = "S12">);</span></p></div></div><p class = "S13"><span class = "S2">We follow the same previous steps to plot the water saturation in our computational domain. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% plot brine saturation on the right side of current figure</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,2);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">PlotCellData(G, S(:), </span><span class = "S18">'EdgeColor'</span><span class = "S17">,</span><span class = "S18">'none'</span><span class = "S12">);</span></p></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S17">axis </span><span class = "S18">tight equal</span><span class = "S17">; colorbar, title(</span><span class = "S18">'Brine saturation [-]'</span><span class = "S12">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 658px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div></div><p class = "S0"></p><div class = 'SectionBlock containment'><p class = "S3"><span class = "S24">Post-processing of the numerical solution - contouring mode</span></p><p class = "S3"><span class = "S2">Presenting the results in mapping mode is most natural when working with cell-centered discretization methods. It has also the advantage to appreciate the dispersion of the saturation plume as seen in the right plot of the above figure. But, it has other limitations. For instance, the gradient of the plotted variable cannot be easily appreciated in some areas, the presentations graphics quality may deteriorate when saving the figure into raster file formats, and this presentation mode is not suitable for explanation to non-specialists such as managers and decision makers. </span></p><p class = "S3"><span class = "S2">Another presentation mode which is quite popular in all fields of geosciences in </span><span class = "S16">contouring</span><span class = "S2">. The target variable iso-values are traced by a general purpose algorithm to contour the data. Such presentation mode is mainly supported by basic MATLAB routines such as </span><a href = 'https://fr.mathworks.com/help/matlab/ref/contour.html'><span class = "S15">contour</span></a><span class = "S2">, </span><a href = 'https://fr.mathworks.com/help/matlab/ref/contourc.html'><span class = "S15">contourc</span></a><span class = "S2">, </span><a href = 'https://fr.mathworks.com/help/matlab/ref/contourf.html'><span class = "S15">contourf</span></a><span class = "S2"> and </span><a href = 'https://fr.mathworks.com/help/matlab/ref/contour3.html'><span class = "S15">contour3</span></a><span class = "S2">. We will use </span><span class = "S14">PlotContourData</span><span class = "S2"> function which is a simple wrapper to </span><span class = "S14">contourf</span><span class = "S2"> or </span><span class = "S14">contour</span><span class = "S2"> to plot filled or simple contour plots, respectively. This function arguments are similar to those of </span><span class = "S14">PlotCellData</span><span class = "S2"> function except for the optional arguments are property-value pairs which are compatible with those of the </span><span class = "S14">contourf</span><span class = "S2"> or </span><span class = "S14">contour</span><span class = "S2"> commands. </span></p><p class = "S3"><span class = "S2">So, make a figure made from two contour plots of water pressure and saturation. We plot pressure contour lines between 20 and 140 bars with an interval of 20 bars. We select saturation contour lines between 0 and 0.7 with 5% interval. </span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">figure;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">opts.fill = </span><span class = "S18">'off'</span><span class = "S17">;   </span><span class = "S11">% change this option to 'on' for contour filled plots</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% plot water pressure on left side </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,1);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">PlotContourData(G, P(:)/bar,opts,</span><span class = "S18">'LevelList'</span><span class = "S12">,20:20:140);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">colormap(jet(16));</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">axis </span><span class = "S18">tight equal</span><span class = "S17">; colorbar, title(</span><span class = "S18">'Pressure [bar)]'</span><span class = "S12">);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S11">% plot water saturation on right side </span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">subplot(1,2,2);</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">PlotContourData(G, S(:),opts,</span><span class = "S18">'LevelList'</span><span class = "S12">,0:0.05:0.7);</span></p></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S17">axis </span><span class = "S18">tight equal</span><span class = "S17">; colorbar, title(</span><span class = "S18">'Brine saturation [-]'</span><span class = "S12">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 658px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S13"><span class = "S2">To better understand the spatial distribution of brine saturation we also plot its probability distribution reflecting the saturation volumetrics inside the reservoir. To do that we call the </span><a href = 'https://fr.mathworks.com/help/matlab/ref/histogram.html'><span class = "S15">histogram</span></a><span class = "S2"> command for </span><span style="vertical-align:-3"><img src="" width="10" height="18" /></span><span class = "S2"> variable with 7 bins as a second argument and '</span><span class = "S14">probability</span><span class = "S2">' value for the '</span><span class = "S14">Normalization</span><span class = "S2">' option. This means that the height of each bar is the relative number of cell-centered saturation items (number of saturation values in bin / total number of grid cells). The sum of the bar heights is obviously equal to 1.</span></p><div class = 'LineNodeBlock contiguous'><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S12">figure;</span></p></div><div class = 'inlineWrapper'><p class = "S10 lineNode"><span class = "S17">histogram(S,7,</span><span class = "S18">'Normalization'</span><span class = "S17">,</span><span class = "S18">'probability'</span><span class = "S12">);</span></p></div><div class = 'inlineWrapper outputs'><p class = "S10 lineNode"><span class = "S17">title(</span><span class = "S18">'Saturation probability distribution'</span><span class = "S12">);</span></p><div class="outputParagraph"><div class="inlineElement embeddedOutputsFigure" style="max-height: 800px; width: 658px;"><div class="figureElement"><img class="figureImage" draggable="false" src=""></div></div></div></div></div><p class = "S13"><span class = "S2">Around 35% of the reservoir volume is filled with water saturations between 0.4 and 0.5 while around 75% of the reservoir volume is filled with water saturations between 0.3 and 0.6 which are clearly the dominant modes of saturations after 50 days of oil drainage by water injection. </span></p></div><p class = "S0"></p><div class = 'SectionBlock containment'><h2 class = "S8"><span class = "S2">Conclusions</span></h2><p class = "S3"><span class = "S2">In merely a little bit more than 100 lines of MATLAB commands (including blank and comment lines) we were able in this tutorial to build a script performing the following tasks: </span></p><ul class = "S5"><li class = "S6"><span class = "S0">Read petrophysical data from disk files, plot, and analyse them.</span></li><li class = "S6"><span class = "S0">Setup auxilliary simulation data.</span></li><li class = "S6"><span class = "S0">Run an immiscible two-phase flow drainage simultion involving an implicit sequential simulation procedure, analyse and plot the convergence behavior of the adaptive time stepping algorithm. </span></li><li class = "S6"><span class = "S0">Post-process the numerical solution using mapping and contouring modes, then analyse the spatial distribution of water saturation after 50 days of oil drainage. </span></li></ul><p class = "S3"><span class = "S2">All these steps take less than 1 minute CPU time in a commodity PC. With a scripting technology the time traditionally taken in file and data transfer between pre-processing, simulation, and post-processing tools as in existing workflows is cutted to its minimum. Hence, the modeller can spend more time to gain in problem understanding rather than loosing a significant project timeframe in useless data manipulation tasks.</span></p></div></div>
<!-- 
##### SOURCE BEGIN #####
%% Immisicble two-phase flow drainage in a subsurface reservoir - example 1 
% In this tutorial two-dimensional immiscible two-phase flow drainage of oil 
% by water injection is simulated. Permeability and porosity distributions are 
% extracted from the first layer of the full three-dimensional geological model 
% of <http://www.spe.org/web/csp/datasets/set02.htm SPE10 comparative solution 
% project>. 
% 
% This tutorial involves four steps:
% 
% # First, permeabiiity and porosity data are loaded from disk files, analysed, 
% and visualised. Next, we verify that these variables are highly correlated as 
% expected from the data.
% # Setup auxilliary simulation data, such as wells positions, their flow rates, 
% fluid phases properties, initial pressure and saturation conditions, the time 
% stepping parameters, and the nonlinear saturation solver options.
% # Run the time simulation loop involving a sequential implicit solution procedure 
% between pressure and saturation distributions. Next, we analyse and plot the 
% convergence behavior and the performance of the nonlinear saturation solver.
% # Post-process of the numerical results by plotting the spatial distributions 
% of simulated pressure and saturation at the simulation end. Different presentation 
% graphics styles such as _mapping_ and _contouring_ modes are explored. 
%% Load and plot the petrophysical (permeability and porosity) data
% We considere a simple one layer rectangular domain. In this case we explicitly 
% set the domain lengths along each space dimension as follows

% Domain size along X, Y & Z directions
Dx = 365.76; Dy = 670.56; Dz = 0.6096;
%% 
% Next, we set the number of cells along each space direction and the total 
% number of grid cells. Next, we add those variables into the |Grid| structure

Nx = 60;
Ny = 220;
Nz = 1; 
N  = Nx*Ny*Nz;
Grid.Nx = Nx; Grid.Ny = Ny; Grid.Nz = Nz; Grid.N = N;
%% 
% To fullfill the need of internal compuattional tasks we add the uniform 
% grid spacings, |hx|, |hy|, and |hz|, along the respective directions to the 
% |Grid| structure

Grid.hx = (Dx/Nx); Grid.hy = (Dy/Ny); Grid.hz = (Dz/Nz);
%% 
% Additionally, we need to setup other few variables for the |Grid| structure 
% such as, |V| (the uniform cells volume), and |compr| (the total compressibility). 
% Notice that the compressibility data is a 3D shape array whose size conforms 
% to the grid dimensions. The MATLAB command |<https://fr.mathworks.com/help/matlab/ref/ones.html 
% ones>(Nx,Ny,Nz)| constructs such an array with a unit values for all cells, 
% and thus by multiplying it with a constant values we obtain a uniform distribution 
% of the compressibility data over the computational domain. 

Grid.V     = (Dx/Nx)*(Dy/Ny)*(Dz/Nz);
Grid.compr = 4.4e-10.*ones(Nx,Ny,Nz);
%% 
% Now, we load porosity data from a previously saved mat file '_spe_layer1_phi.mat_' 
% which is included in the data subfolder of the included files. porosity data 
% is stored in _phi_layer_ variable (because we can store more than one array 
% in a mat file).

load 'spe_layer1_phi.mat' phi_layer;
%% 
% Now, we will check if the loaded data is within acceptable bounds. A situation 
% which occur many often in practice. First, look if there exists porosity values 
% below or close to 0+. This can be simply done with the following command:

length(find(phi_layer<eps))
%% 
% indicating that the porosity is below the smallest floating point precision 
% number |<https://fr.mathworks.com/help/matlab/ref/eps.html eps>| (effective 
% zero on the target computer) at 68 cells. To correct for this misleading data 
% we set these values to a prescribed minimum porosity threshold (1e-3 herein). 
% Otherwise, the saturation solver will fail since dividing by zero porosity will 
% result in an overflow.

phi_layer = max (phi_layer(:), 1e-3);
%% 
% The |Grid| structure is promoted with a new variable |por| conforming 
% to the |Grid| shape. 

Grid.por = reshape(phi_layer',Nx,Ny,Nz); 
%% 
% Finally, permeability data is loaded from the mat file '_spe_layer1_perm.mat_' 
% and variable _perm_layer_. 

load 'spe_layer1_perm.mat' perm_layer;
%% 
% Tha data units is converted into SI unit of the permeability (m^2) because 
% it is given in Dary units in the SPE10 dataset. Again, the |Grid| structure 
% is promoted with a new variable |K| conforming to the |Grid| shape but with 
% 3 components since the permeability data assumes an anisotropic diagonal tensor. 

Grid.K = perm_layer*9.869232667160130e-13;
Grid.K = reshape(Grid.K',3,Nx,Ny,Nz); 
%% 
% Now, we're ready the plot the porosity and permeability data to see how 
% they look like. First, we will plot the porosity data in linear scale. Next, 
% we plot permeability along X direction in log10 logarithmic scale to depict 
% its extent over several orders of magnitude. 
% 
% To do this, a cartesian grid |G| data structure is setup from the domain 
% extent ans spacings. This will help to do the plotting commands more efficiently.   

G = cartGrid(0:Dx/Nx:Dx, 0:Dy/Ny:Dy);
%% 
% Start a new MATLAB figure and specify the colormap corrsponding the color 
% palette to use and the number of colors in this palette.

figure;
colormap(jet(16));
%% 
% Plot the porosity data into a |<https://fr.mathworks.com/help/matlab/ref/subplot.html 
% subplot>| in the same window but at le left hand side. The |PlotCellData| function 
% plots a map of the variable |Grid.por| using information in the structure variable 
% |G|. The '|EdgeColor|' option is set to '|none|' to hide the display o the gridlines 
% of the underlying grid. Since the grid is so dense in our modl, letting this 
% option active the resulting plot will be obsured. 

subplot(1,2,1);
PlotCellData(G, Grid.por(:), 'EdgeColor','none');
%% 
% next, this plot x and y axis are set as dependent and we add a color bar 
% and a title. 

axis tight equal; colorbar, title('Porosity map [-]');
%% 
% The same sequence of commands are repeated for the right plot to obtain 
% a map of Kx in log10 scale. 

subplot(1,2,2);
PlotCellData(G, reshape(log10(Grid.K(1,:,:,:)),N,1), 'EdgeColor','none');
axis tight equal; colorbar, title('log_{10}(K_x) permeability map [m^2]');
%% 
% Clear phi_layer variable from memory space since it's not needed further. 
% This is more relevant for large scale applications or on computers with very 
% modest ressources. 

clear phi_layer;
%% 
% Now, we check if the log10-permeability and porosity fields ara correlated 
% somehow. First, let's start with the first component of the diagonal permability 
% tensor, Kx:

corrcoef(reshape(log10(Grid.K(1,:,:,:)),N,1),reshape(Grid.por(:,:,:),N,1))
%% 
% Off-diagonals of the correlation matrix are high positive values indicating 
% a good positive correlation between the two variables. Let's now do the same 
% thing for the y and z components of the permeability to conclude on their correlation 
% with the porosity. 

corrcoef(reshape(log10(Grid.K(2,:,:,:)),N,1),reshape(Grid.por(:,:,:),N,1))
corrcoef(reshape(log10(Grid.K(3,:,:,:)),N,1),reshape(Grid.por(:,:,:),N,1))
%% Setup auxilliary simulation data
% We further add a five-spot well pattern composed from 1 injector and 4 producers. 
% Production wells are placed at the domain corners while the injection well is 
% located at the center. A quarter of the volumetric injection flow rate, |Qi|, 
% is produced at each of the four production wells. 

% define some useful time units conversion constants 
second = 1; 
minute = 60*second;
hour   = 60*minute; 
day    = 24*hour; 
week   = 7*day;
month  = 30*day; 
year   = 365*day;

% Injection/production flow rates
Qi         = 25/hour;                  % total flow rate (converted from m^3/h to m^3/s)
prod_wells = [1,Nx,N-Nx+1,N];          % grid indices of production wells
Qw         = zeros(N,1);               % initialize production/injection flow rates array
Qw(prod_wells) = -Qi/4*ones(1,4);      % flow rates at production wells 
inj_well   = (Nx*(Ny-1))/2;            % grid index of injection well
Qw(inj_well)   = +Qi;                  % injection flow rate
%% 
% Notice that (by convention), flow rates have to be positive for injection 
% wells and negative for production wells. |Qw| is a distributed 1xN array which 
% holds distributed flow rates. It has zero values everywhere except at injection/production 
% cells. 
% 
% Next, we initialize two fluid objects |oil| and |brine| in which fluid 
% properties data corresponding to a more dense and viscous _Oil_ phase and less 
% viscous and dense _Saltwater_ phase, respectively. 

% fluid properties units: density [Kg/m^3], viscosity [Pa.s], residual saturation []
oil   = Fluid('Oil',      [1200, 3e-3,  0.3]);   % resident fluid
brine = Fluid('Saltwater',[1030, 5e-4, 1e-2]);   % injected fluid
%% 
% Initialize fluid global pressure in the computational domain to a constant 
% value of 100 bars and the initial brine saturation to its residual value (i.e. 
% |brine.Sr|)

% define some useful pressure units conversion constants
Pascal = 1;                                      % Pascal SI unit 
bar    = 1e+5*Pascal;                            % bar
psi    = 7.042253521126761e+03*Pascal;           % UK pressure unit = lbf/in^2
atm    = 1.013249965828145e+05*Pascal;           % atmosphere eqv unit 

P = 100*bar*ones(Nx,Ny,Nz);
S = brine.Sr.*ones(N,1);
%% 
% Define time stepping parameters including total simulation time, |T|, 
% the number of time steps, |nt|, and the global time step, |dt.|

T   = 50*day;                      % 50 days
nt  = 25;                          % number of time steps
dt  = T/nt;                        % time step
%% 
% For the nonlinear solvers we need to pass an optional structure with options 
% holding the absolute tolerance, tol, and the maximum number of allowed Newton-Raphson 
% iterations, maxiter. The latter corresponds to the inner iteration count since 
% the most outer loop corresponds to time stepping. The last option, min_dt, is 
% the minimum allowed time step size for the saturation solver. This is needed 
% because, as explained later in more detail, we use an adaptive time stepping 
% algorithm which cancels out if the time step size becomes smaller than this 
% prescribed threshold. 

opt.tol     = 1e-5;
opt.maxiter = 100;
opt.min_dt  = 30*second;
%% Run the time simulation loop - Sequential implicit solution
% Now, we will build the time loop involving a sequential solution of the governing 
% global pressure equation and the nonlinear non-wetting phase saturation equation 
% formulated in fractional form. The pressure and saturation solvers are both 
% implicit implying an unconditional stability regardlless of the magnitude of 
% the selected time step. 
% 
% The pressure equation, neglecting the gravitational and capillary pressure 
% terms is given as:
% 
% $$\phi c_t \frac{\partial P}{\partial t} = \nabla.[K \lambda_t \nabla p] 
% + \frac{q_w}{\rho_w} +  \frac{q_o}{\rho_o}$$ 
% 
% where $c_t$ is the total compressibility, $\phi$ is the rock porosity, 
% $P$ is the fluid pressure, $\lambda_t = \lambda_o + \lambda_w$ is the total 
% mobility given as a sum of oil and water mobilities, respectively. A fluid phase 
% mobility $\lambda_l$ is defined as the ratio of its relative permeability $k_{r,l}$ 
% to its viscosity $\mu_l$. $q_w$ is the water injection/withdawal mass flow rate 
% [M/T] and $\rho_w$ is the density of the water phase [M/L^3]. Threfore $\frac{q_w}{\rho_w}$ 
% is the volumetric flow rate.
% 
% The saturation equation is a nonlinear diffusion-convection like equation 
% for the water phase given as:
% 
% $$\phi \frac{\partial S_w}{\partial t} = - \nabla . (f_w \textbf{u}) + 
% \frac{q_w}{\rho_w}$$
% 
% where $S_w$ is the water phase saturation, $f_w(S_w)$ is the saturation-dependent 
% fractional flow function for the water phase, $\textbf{u} = \textbf{u_{w} + 
% \textbf{u_o}$ is the total velocity (i.e. the sum of fluid velocities of each 
% phase). 
% 
% Herein, the pressure solver is linear but requires successive updating 
% of the total mobility of both fluid phases which is the task of |TwoPhasePressure 
% |function. The latter wraps the _Two-Point Flux Approximation solver_ (|TPFA| 
% function) such that $\lambda_t (S_w)$ is iteratively updated from spatial distribution 
% of water saturation calculated in the previous time step. 
% 
% The saturation solver uses an heuristic adaptive time steeping procedure 
% combined with an inner Newton-Raphson (NR) iteration loop. During the inner 
% NR iterations, the sparse system Jacobian matrix $\textbf{J} = [\frac{dF(S_w)}{dS_w}]$ 
% is evaluated analytically to accelerate the algorithm convergence. If the NR 
% iteration fails the current time step is subdivided (i.e. the number of local 
% time steps is multiplied by 2) and the whole time step is repeated again from 
% the beginning. The latter mechanism is called back-stepping. The local NR failure 
% is due to either (i) an insuffisant or small number of allowed maximal iterations 
% to achieve the required accuracy. Hence, when setting the NR options the user 
% should keep a good balance between the |maxiter| and |tol| options; (ii) second, 
% to a high or very difficult level of prescribed accuracy which cannot be achived 
% for badly preconditioned problems. So, the absolute tolerance should be kept 
% to reasonable values (in practice values between 1e-3 and 1e-4 are suggested); 
% (iii) and finally, when the selected global time step results in large local 
% pressure jumps, the same time step is generally inappropriate for the saturation 
% solver implying a smaller time step advancement of the saturation solver to 
% absorp the sudden numerical shock caused by the large pressure jump. 
% 
% In general, this adaptive time stepping strategy is effective in practice. 
% Very often the algorithm starts with smaller saturation time steps at the begining 
% of the simulation when initial conditions are far from the expected solution. 
% This sub-time step will gradually increase in the next time levels if there 
% is no abrupt changes in the boundary conditions (i.e. wells pumping or production 
% schedules) or materials and fluids properties. 
% 
% Let's run this two-phase flow oil drainage simulation

tic;              % start timing the time marching loop
for t=1:nt
   fprintf('\n'); 
   fprintf('Solving two-phase flow problem. Time = %f days\n', t*dt/day);
   
   % call TPFA flow solver adapted to two-phase flow 
   [P,V] = TwoPhasePressure(Grid,S,brine,oil,Qw,P,dt); 
   
   % solve for brine saturation 
   [S,out(t)] = ImplicitSaturation(Grid,S,brine,oil,V,Qw,dt,opt);
   
   % quit any time step if the saturation solver fails 
   if out(t).convergence==false, break; end
end
toc;
%% 
% From the convergence output we notice that this sequential solution algorithm 
% went very smoothly. There was a difficulty to reach convergence in the first 
% time step because it seems that 2 days time step was very large for the saturation 
% solver. Hence, the initial time step was automatically subdivided 11 times resulting 
% to an effective saturation time step close to 84 seconds as shown below from 
% the fields of the returned output structure |out| for the first time period. 
% The |dt| field returns the selected time step by the saturation solver during 
% this period. 

out(1)
%% 
% Now, we plot a summary of the saturation solver performance during the 
% simulation. In the left side of the figure we plot the evolution of the internally 
% selected time step versus time, while in the right figure the evolution of the 
% total number of inner Newton-Raphson iterations versus time is plotted. Both 
% figure are best depicted using a |<https://fr.mathworks.com/help/matlab/ref/semilogy.html 
% semilogy>| axis. 

figure; 

subplot(1,2,1); 
for t=1:nt, dt_history(t) = out(t).dt; end 
semilogy((1:nt)*dt/day, dt_history/day, '-ro'); 
title('Evolution of the saturation time-step'); xlabel('Time [days]'); ylabel('Time step [days]'); grid on;

subplot(1,2,2); 
for t=1:nt, nr_history(t) = out(t).nr_iterations; end 
semilogy((1:nt)*dt/day, nr_history, '-b*'); 
title('Newton-Raphson iterations'); xlabel('Time [days]'); ylabel('Total NR iterations'); grid on;
%% 
% We observe thet the internal saturation time step gradually increases 
% to 1.5 hour, 6 hours, 12 hours and finally to 2 days starting from 22 days after 
% the simulation start. So, afterwards the pressure and saturation time steps 
% becomes equal meaning that the saturation extent of the plume becomes significant 
% enough to amortize the pressure jumps immediately. 
% 
% *Post-processing of the numerical solution - mapping mode*
% 
% Since the main simulated variables are pressure $P$ and water saturation 
% $S_w$ we will plot them directly using the supplied functions which make use 
% of MATLAB's excellent built-in graphics capabilities without writing the results 
% into intermediate text or binary files in a required format and loading them 
% back from an external third-party post-processing software. So, with very few 
% commands we get directly the results ready for quality publication freeing the 
% modeller from all intermediate tasks and hence speeding the data-processing-simulation-post-processing 
% workflow significantly. 
% 
% We will use here the routine |PlotCellData| to produce the figures in _mapping_ 
% mode where each cell in the grid is plotted with a constant color corresponding 
% to a linearly interpolated value between the upper and lower colors in the used 
% color palette. This palette fits the range between the maximal and minimal values 
% detected in the data array.  
% 
% The |<https://fr.mathworks.com/help/matlab/ref/figure.html figure>| command 
% without any argument opens a new MATLAB figure (or axis in MATLAB's terminology) 
% so that the previously opened figures are preserved and all following drawing 
% commands will be placed in this new figure. The |<https://fr.mathworks.com/help/matlab/ref/colormap.html 
% colormap>| command is optionally called before the drawing commands to choose 
% among the available color palettes, here we choose the |<https://fr.mathworks.com/help/matlab/ref/jet.html 
% jet>| palette with 256 colors distributed between blue and red. 

figure; 
colormap(jet(256));
%% 
% Next, we plot the pressure after conversion from |Pascal| to |bar| units. 
% We call the |<https://fr.mathworks.com/help/matlab/ref/subplot.html subplot>| 
% command to specify that we want to make a plot matrix with 1 row and 2 columns 
% in the target figure. The last argument in the subplot command is the active 
% subplot to which the next plotting commands are directed. The |PlotCellData| 
% function takes additional property-value pairs (à la MATLAB). These options 
% are exactly those supported by the low level |<https://fr.mathworks.com/help/matlab/ref/patch.html 
% patch>| command to control the internal details of the plot. To study all these 
% options type |help patch| command or go to |<https://fr.mathworks.com/help/matlab/ref/patch.html 
% patch>| window in the online help system.  Here the '|EdgeColor|' property specifies 
% the color of the underlying grid edges, by setting this property to '|none|' 
% they will be hidden in the current plot. In this example the grid is very dense. 
% So, drawing all grid edges will obscure the reulting plots but in other applications 
% where the grid is coarse this might not be necessary.

% plot fluid pressure on left side of current figure
subplot(1,2,1);
PlotCellData(G, P(:)/bar, 'EdgeColor','none');
%% 
% By default MATLAB plots fit to the rectangular window inside which they're 
% displayed. A good idea here is to call the |<https://fr.mathworks.com/help/matlab/ref/axis.html 
% axis>| command to make x and y axis scales dependent using |equal |option. Likewise, 
% |tight| wil readjust the plot to fit at best to actual window after axis scaling. 
% Finally, we decorate the plot with a |<https://fr.mathworks.com/help/matlab/ref/colorbar.html 
% colorbar>| and a custom |<https://fr.mathworks.com/help/matlab/ref/title.html 
% title>|. 

axis tight equal; colorbar; title('Pressure [bar]');
%% 
% We follow the same previous steps to plot the water saturation in our 
% computational domain. 

% plot brine saturation on the right side of current figure
subplot(1,2,2);
PlotCellData(G, S(:), 'EdgeColor','none');
axis tight equal; colorbar, title('Brine saturation [-]');
%% 
% *Post-processing of the numerical solution - contouring mode*
% 
% Presenting the results in mapping mode is most natural when working with 
% cell-centered discretization methods. It has also the advantage to appreciate 
% the dispersion of the saturation plume as seen in the right plot of the above 
% figure. But, it has other limitations. For instance, the gradient of the plotted 
% variable cannot be easily appreciated in some areas, the presentations graphics 
% quality may deteriorate when saving the figure into raster file formats, and 
% this presentation mode is not suitable for explanation to non-specialists such 
% as managers and decision makers. 
% 
% Another presentation mode which is quite popular in all fields of geosciences 
% in _contouring_. The target variable iso-values are traced by a general purpose 
% algorithm to contour the data. Such presentation mode is mainly supported by 
% basic MATLAB routines such as |<https://fr.mathworks.com/help/matlab/ref/contour.html 
% contour>|, |<https://fr.mathworks.com/help/matlab/ref/contourc.html contourc>|, 
% |<https://fr.mathworks.com/help/matlab/ref/contourf.html contourf>| and |<https://fr.mathworks.com/help/matlab/ref/contour3.html 
% contour3>|. We will use |PlotContourData| function which is a simple wrapper 
% to |contourf| or |contour| to plot filled or simple contour plots, respectively. 
% This function arguments are similar to those of |PlotCellData| function except 
% for the optional arguments are property-value pairs which are compatible with 
% those of the |contourf| or |contour| commands. 
% 
% So, make a figure made from two contour plots of water pressure and saturation. 
% We plot pressure contour lines between 20 and 140 bars with an interval of 20 
% bars. We select saturation contour lines between 0 and 0.7 with 5% interval. 

figure;
opts.fill = 'off';   % change this option to 'on' for contour filled plots
% plot water pressure on left side 
subplot(1,2,1);
PlotContourData(G, P(:)/bar,opts,'LevelList',20:20:140);
colormap(jet(16));
axis tight equal; colorbar, title('Pressure [bar)]');
% plot water saturation on right side 
subplot(1,2,2);
PlotContourData(G, S(:),opts,'LevelList',0:0.05:0.7);
axis tight equal; colorbar, title('Brine saturation [-]');
%% 
% To better understand the spatial distribution of brine saturation we also 
% plot its probability distribution reflecting the saturation volumetrics inside 
% the reservoir. To do that we call the |<https://fr.mathworks.com/help/matlab/ref/histogram.html 
% histogram>| command for $S$ variable with 7 bins as a second argument and '|probability|' 
% value for the '|Normalization|' option. This means that the height of each bar 
% is the relative number of cell-centered saturation items (number of saturation 
% values in bin / total number of grid cells). The sum of the bar heights is obviously 
% equal to 1.

figure;
histogram(S,7,'Normalization','probability');
title('Saturation probability distribution');
%% 
% Around 35% of the reservoir volume is filled with water saturations between 
% 0.4 and 0.5 while around 75% of the reservoir volume is filled with water saturations 
% between 0.3 and 0.6 which are clearly the dominant modes of saturations after 
% 50 days of oil drainage by water injection. 
%% Conclusions
% In merely a little bit more than 100 lines of MATLAB commands (including blank 
% and comment lines) we were able in this tutorial to build a script performing 
% the following tasks: 
% 
% * Read petrophysical data from disk files, plot, and analyse them.
% * Setup auxilliary simulation data.
% * Run an immiscible two-phase flow drainage simultion involving an implicit 
% sequential simulation procedure, analyse and plot the convergence behavior of 
% the adaptive time stepping algorithm. 
% * Post-process the numerical solution using mapping and contouring modes, 
% then analyse the spatial distribution of water saturation after 50 days of oil 
% drainage. 
% 
% All these steps take less than 1 minute CPU time in a commodity PC. With 
% a scripting technology the time traditionally taken in file and data transfer 
% between pre-processing, simulation, and post-processing tools as in existing 
% workflows is cutted to its minimum. Hence, the modeller can spend more time 
% to gain in problem understanding rather than loosing a significant project timeframe 
% in useless data manipulation tasks.
##### SOURCE END #####
--></body></html>